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ABSTRACT 


A numerical study of thermal convection driven by uniform volumetric 
energy sources has been conducted for a horizontal fluid layer bounded 
from above by a rigid, isothermal stirface and from below by a rigid, zero- 
heat-flux surface. The side walls of the fluid domain are assumed to be 
rigid and perfectly instilating. The computations are formally restricted 
to two-dimensional laminar convection but have been carried out for a 
range of Rayleigh numbers which spans the regimes of laminar and turbu- 
lent flow . 

The results of the computations consist of strea m line and isotherm 
patterns, horizontally averaged temperature distributions, and horizon- 
tally averaged Nusselt numbers at the upper surface. Flow and temperature 
fields do not exhibit a steady state, but horizontally averaged Nusselt 
numbers reach limiting, quasi-steady values for all Rayleigh numbers con- 
sidered. Correlations of the Nusselt number in terms of the Rayleigh and 
Prandtl nmbers have been determined in the following forms : 

Nui = 0.477 Ra°*^^° Pr0.04Cf7 

5 X 10^ g Ra S 5 X 10® 

0.05 g Pr g 20, 


and 


,Nui = 0.420 
5 X 10® g Ra S 5 X 10® 

Pr = 6 . 5 . 

The latter correlation is in excellent agreement with existing experi- 
ments . Horizontally averaged tanperature distributions are in good agree- 
ment with time-averaged temperatvire measurements in the layer. 

■ ^^Streamline and isotherm patterns for Pr = 6.5 indicate that the 
\'Gohve"cttye'prbce'ss is dominated by up-flows of warm fluid over broad re- 
gions "of the layer'. Faster down-flows are confined to narrow regions, 
and dovm-flows are always observed on the side walls of the fl\iid domain. 
At ’low Raylpigh numbers, up-flows occ\ir in the center of the layer, and 
two regions of recirculating flow are observed. At higher Rayleigh num- 
bers, down-flow is observed at the center of the layer with four regions 
of recirculating flow occupying the layer. For Prandtl numbers other 
than 6 . 5 , the streamline and isotherm patterns are altered, and the fea- 
tures of these flows are compared with the Pr = 6.5 cases for selected 
Rayleigh numbers. 
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SECTION I - INTRODUCTION 


The ptirpose of this study is to compute via finite-difference meth- 
ods the temperature and flow fields of two-dimensional, thermal convection 
which is driven by uniform voltimetric energy sources in a horizontal fluid 
layer. It is assumed that fluid satisfies the Boussinesq equation of 
state and is a normal fluid in the sense that p > 0. The fltoid layer is 
taken to be bo-unded from above and below by rigid surfaces which are of 
constant temperature and zero heat flux, respectively. Vertical bound- 
aries of the layer are taken to be perfectly insxalating, and the layer 
aspect ratio (i.e., depth divided by horizontal extent) is varied from 
0.125 to 1.0. This variation of the layer aspect ratio is, perhaps, su- 
perfluous because the aspect does not appear explicitly in the mathemat- 
ical formulation of the problem; however, such calculations are performed, 
in part, as a check on the computational algorithm and for the sake of 
completeness. Computations are also done for the case of a layer bounded 
from above and below by two rigid, isothermal surfaces owing to the exist- 
ence of both experimental [ 7 ] and theoretical [12] studies with these 
boundary conditions. 

Additional objectives of the present study are to investigate the 
influence of the Prandtl number on the flow and temperature fields in the 
layer and to compute the heat fl\ix at the upper surface of the layer. 

Owing to both stability and convergence criteria and limitations on com- 
putational time, the Prandtl number range considered in the present study 
is 0.05 S Pr S 20. Heat transfer results are presented as correlations 
in terms of the Nusselt nvimber at the upper stirface and the Rayleigh and 
Prandtl numbers in a form suitable for engineering applications. 

From the known theoretical works on thermal convection, it appears 
that exact solutions for the temperature and flow fields can be obtained 
only for a relatively small range of Rayleigh numbers via either approxi- 
mate or exact analytical methods [I 5 ]. On the other hand, the applica- 
tion of finite-difference methods to approximate the system of partial 
differential equations governing the convective process can yield solu- 
tions for a larger range of Rayleigh numbers. Numerical work up to the 
present has treated thermal convection with uniform volumetric energy 
sources in a horizontal fluid layer, with thermal and hydrodynamic bound- 
'a.vy conditions of interest here, only for low Rayleigh number laminar con- 
vection, e.g., the numerical work by Thirlby [I 6 ] for I .5 x 10® i 
Ra ^ 2.6 X lO'^. Therefore, our goal is to obtain numerical solutions for 
the flow and temperature fields for Rayleigh numbers larger than those of 
previous studies. 

The computational approach taken in the present study is to cast the 
governing time-dependent nonlinear partial differential equations in 
finite-difference form and to seek quasi- steady solutions in terms of a 
convergence criterion placed on the temperature distribution within the 
layer. This approach is followed in recognition of studies of Kulacki 
and Goldstein [7], Maylnger, Jahn, Reineke and Steinberner [12], and 
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Goldstein, Chu and Kulacki [6] which showed both experimentally . and the- 
oretically that in thermal convection with internal energy sources in a 
layer with two constant temperature boundaries, horizontally averaged 
two-dimensional temperature and flow fields exhibit no truly steady pat- 
tern despite the fact the average heat fluxes at the layer boundaries 
reach steady values at a given Rayleigh number. It was expected that the 
same behavior would be foimd in the present study because the convective 
and momentum and energy transport processes within the layer are not much 
affected by the substitution of a zero heat flux surface at the lower 
boundary. 
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SECTION II - GOVERNING EQUATIONS 


The temperature and flow fields of thermal convection are governed 
mathematically by a coupled system of nonlinear partial differential 
equations derived from the conservation requirements for mass, momentum 
and energy [2]. In these equations, the Boussinesq approximation [13] 
is used, i.e., all of the thermophysical properties are assumed to be 
constant except the density in the body force term which is related to 
temperature by 


p = PrCl - p(T - Tj.)] , 


( 1 ) 


where pj. is the density corresponding to the temperature Tp. The conser- 
vation equations for mass, momentvim, and energy are, 


V • u = 0 , 

Su V 

^ + (u • V) u = — ^ + g[l - P(T - Tp)l + vV^u , and 
at p 

+ (u . V) T = C^T + ^ . 
dt - pC 

By introducing the variables 

0 = T - Tp, Tp = Ti , 

P = p + pgy , and 

S = H/pC , 

the momentum and energy equations become 

VP 

^ + (u • V) u + Pg0J^ + vV^u 


( 2 ) 

(3) 

( 4 ) 

(5) 

( 6 ) 

(7) 

(8) 


+ (u . V) 0 = q!V^ 0 + S . 
dt " ■ 


(9) 


The conservation equations may be expressed in dimensionless form by 
introducing the following variables: 
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V = LV , 
ta 


t* = 




u' = 


p« = 


u 

a/L ’ 


0 ' 


pva 

0 


, and 


(!)■ 


Equations (2), (8), and (9) can then be written 

V • u = 0 , 


^ + (u • V) = -VP + Ra 0j_ + V^u , and 


3t 


+ (u • V) 0 = v^e + 2 , 


( 10 ) 


( 11 ) 

( 12 ) 

(13) 


where the primes have been dropped in Eqs.. (3JL-13). 

By taking the curl of the momentum equation, Eq. (12), the pressure 
term is eliminated under consideration of the identity Vx(VP) = 0, and 
Eq. (12) can be written 

+ (u . V) «) = -Ra k + , (14) 

where the vorticity, cj = Vxu, has been introduced. 


In the case of two-dimensional processes, the vorticity has only one 
component and therefore it can be treated as a scalar quantity. Therefore, 
the energy and momentum equations can be written as 


80 8(u0) 8(v0) _ ^0 ^0 p 

St Sx Sy Sx:^ Sy^ ’ 

1 /Su S(uoj) . S(vw)\ „ S0 , S^CJ , 

Sx ^ ^ ^ ^ ^ ’ 

where the vorticity is defined as 

Su Sv 


(15) 

(16) 


(17) 
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and the continuity equation, Eq. (ll) , has been combined with Eqs. (I 3 ) 
and (l4) . 


In order to calculate the velocity field, u = (u,v), from the vor- 
ticity field, it is useful to introduce the stream function, \|r, which is 
related to the velocity field by 

u.|,v=-|. ( 18 ) 


From Eqs. ( 17 ) and (I 8 ) , the relation between the vorticity and stream 
functions can be written as 


0) 




(19) 


Boundary and initial conditions on u and 0 are next needed to com- 
plete the mathematical formulation of the problem posed by Eqs. (15-19)* 

It is assumed that fluid is contained in a rectangular domain whose sides 
and bottom are insulated and its top is maintained at constant temperature 
(Fig. 1). Furthermore, all of the boundaries of the domain are assvnned 
to be rigid. Therefore, the thermal and hydrodynamic boundary ; conditions 
associated with Eqs. (15-19) can be form\alated as the following; 


X = 0, X = X/L 

u = V = ijf = 0, S0/Sx = 0 

(20) 

y = 0 

u = V = \|/ = 0, ^ 0 /^ = 0 

(21) 

y = 1 

u=v=)j/ = 0, 0 = 0 

(22) 


The layer is assumed to be at a constant 
t = 0. Thus, the initial conditions are 

u(x,y,o) = v(x,y,0) 

t(x,y,0) = w(x,y,0) 

0(x,y,o) = 0 . 


temperature and motionless at 


= 0. 

= 0 , and V 


(23) 


« 

It should be mentioned that additional boundary conditions are needed 
to obtain the vorticity at the walls. Since the value for the vorticity 
at the rigid walls is difficult to obtain in an explicit form, an approx- 
imate method will be used and this will be discussed in more detail in 
Section 3-2. 
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SECTION III - THE FINITE-DIFFERENCE SOLUTION 


3-1 Fommlation of the Finite-Difference Equations 


The mathematical description of thermal convection with uniform 
volumetric energy soxirces in a horizontal fluid layer by means of coupled, 
nonlinear partial differential equations makes the exact solution for 
such a system intractable. Finite-difference methods are of practical 
value, however, for obtaining an approximate solution to the problem spec- 
ified by Eqs. (15-23). In order to obtain a solution via finite-differ- 
ences, the domain of interest is covered with a grid net, and the flow 
and temperature fields are evaluated only at the grid points at a given 
point in time. In this study, a grid network is obtained by constructing 
a series of equally spaced vertical and horizontal lines parallel to the ' 
X- and y-axes (Fig. l) . The subscripts i and j are used to denote the 
positions of the grid points; i.e., x = (i - 1) Ax, y = (d - l) Ay* 


In order to transform Eqs. (15-23) into finite-difference form, a 
Taylor's series expansion is used to approximate derivatives at a point 
in terms of the function (e.g., temperature or vorticity) at that point 
and its neighboring points . Using a Taylor ' s series expansion for an ar- 
bitrary function f(x,y), the following relations can be written: 


f‘i+l,j = 


fi,4 + («) 


(Ax:)^ Sff 

i,j 2! ^ 




, (Ax)^-^ 

(m-1)! Sx^-^ 


+ R 


m > 


(24) 


and 


^f 


fi-l,J = fi,j - (ax) ^ 


1 , J 


. (Ax)^ ^ 
2! Sx2 




+ (-Ax)”^~^ 


(m-l)! ^3^^' 


-1 


+ R„ 


(25) 




where R„ is the remainder term. From Eqs. (24) and (25), the following 


be 

written: 



Sf 

^i+l,j ■ ^i,J 

+ 0(Ax) , 

( 26 ) 

9x 

i,j Ax 

5f 

^i,.1 - ^i-1, j 

+ 0(Ax) , 

(27) 


i,J Ax 
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Sf 

Sx 


(28) 


I 

8 ? 


^ , and 

i,J 2 ax ^ ' 

(Ax) 


(29) 


1 , J 


Equations (26) and (27) are forward and backward difference approxi- 
mations, respectively, while Eqs. (28) and (29) are central difference 
approximations. The last terms in Eqs. (26-29) represent the truncation 
error and, thus, the acciiracy of the finite-difference approximation. It 
is obvious that the central difference approximations are of higher order 
accuracy than the other representations. However, from the standpoint of 
stability requirements, which will be discussed in Section 3-3 below, the 
representation of the first-order derivatives in Eqs. (l5) and (l6) by 
the central difference formulation is useful only for cases where the 
Rayleigh number is low. In such cases, small grid sizes must be used in 
order to satisfy the stability reqmrement . For example, the computations 
done in this study show that the dimensionless velocities, u and v, reach 
maximum values on the order of 500 for Ra = 5 x 10®, which requires Ax 
and Ay to be less than 0.004. The use of such small grid sizes requires 
a very large storage capacity in the computer, and moreover, a tremendous 
amo\mt of computation time. There is, however, no restriction on the 
grid size when one uses either the backward or the forward difference 
formulation for the nonlinear terms in Eqs. (15) and (l6) . 


In the early work of Barakat and Clark [l], the • nonlinear terms in 
the energy and vorticity equations were approximated with two-point back- 
ward or forward differences according to whether the velocities u and v 
were either positive or negative, respectively. This was termed the 
"first upwind differencing method." Torrance [17] developed a modified 
form of first upwind differencing method in which. 


x-component of velocity u = 


_ ^i,,i ^+l,j 


for example, 

, in the nonlinear term 


the mean 
8u9 
8x 


of 


Eq. (15) was multiplied by or according to whether u was 


positive or negative, respectively. The same method was applied to first 
derivatives with respect to y by replacing the velocity u by v and inter- 
changing i with j. The nonlinear terms in Eq. (16) were similarly treated 
by replacing 0 by co. 


The method developed by Torrance [17] 5 called the "second upwind dif- 
ferencing method"; preserves the conservative and transportive properties 
of the energy and momentum equations. The second upwind difference for- 
mulation of the energy equation, Eq. (15) > and the vorticity equation, 

Eq. (16), both posses the conservative property because the net energy or 
moment-urn transport from the boxindaries of the grid system into its inte- 
rior balances the net increase of the energy or momentum within the sys- 
tem [l4]. In addition, the second upwind differencing method maintains 
the transportive property of the governing equations because any pertvir-r 
batioh in the velocity, energy, or temperat-ure is advected only in the 
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direction of the velocity; i.e., downstream with the flow. All methods 
which use the central difference approximation for the advective terms do 
not possess this property [l4]. Furthermore, Torrance [I7] reported that 
less computation time was needed with the second upwind differencing meth- 
od than with the first upwind differencing method. Also it was shown by 
Torrance that the finite-difference representation of the energy and mo- 
mentum equations using central differences to approximate the advective 
terms (e.g., in the work of Fromm [5l> leads to numerically induced os- 
cillations in the results. Thus, it appears that the use of central dif- 
ferences to approximate the nonlinear terms in Eqs. (15) and (16) can 
lead to computational difficulties. However, the modified scheme for 
forward or backward differences developed by Torrance retains some fea- 
t\ares of central differences, in particular second-order accuracy [l4]. 
Owing to all of these factors, the second upwind differencing method ap- 
pears to be the most suitable approximation for the nonlinear terms in 
Eqs. (15) and (16) and has been adopted in the present work. 

Finite difference methods for solving parabolic partial differential 
equations such as the energy and the vorticity equations, Eqs. (15) and 
(16) respectively, can be classified as either explicit or implicit. The 
time level at which the spatial derivatives are evaluated determines wheth- 
er the scheme is explicit or implicit. If the values of the function at 
the present time, where its values are known at all nodal points, is used 
in Eqs, (26-29), the scheme is termed explicit. The explicit scheme en- 
ables direct computation of the function at all nodal points using a 
simple marching in time procedure. Use of the explicit method, however, 
may require small time increments (i.e., large machine time) in order to 
satisfy stability requirements. To avoid the restriction on the time in- 
crement, implicit methods are usvially recommended. In the implicit meth- . 
od, iterative techniques are generally used for the solution of the 
resulting system of algebraic equations. The Gauss-Seidel iteration pro- 
cedure [19] is suitable for this purpose but may require a large number 
of iterations per time step. On the other hand, increasing the size of 
the time step would increase the number of iterations required to achieve 
any reasonable degree of accuracy. The use of the implicit method from 
the standpoint of realizing a savings of computation time thus may be of 
marginal benefit [17l. Furthermore, the use of the implicit method for 
the solution of the vorticity equation may have a limited advantage over 
the explicit method owing to the lack of a way to explicitly evaluate vor- 
ticity at rigid boundaries. Implicit methods require the use of the vor- 
ticity at boundary nodes at time level n + 1 in order to advance the 
vorticity value at interior nodes from time level n to n + 1. The bound- 
ary value of vorticity at time level n + 1 is not known; and therefore, 
the value of the vorticity at the wall at time level n has to be used to 
approximate that at time level n + 1. Such a linearization of the vor- 
ticity boundary condition requires the use of small time incr^ents so 
that, uP’ at the boTondary will be a good approximation for 

With all of the above factors taken into consideration, the explicit 
formvilation of the finite-difference equations has been used in the pres- 
ent study, with the second upwind differencing method developed by 
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Torrance [ 17 ] to approximate the nonlinear terms in Eqs. (15) and (I6) 
and central differences, Eq. (29), to approximate the diffusion terms. 

The nonlinear spatial derivatives, ^ ^ and , 

ox ay Ox oy 

are approximated with special three-point noncentral differences [l^l] by 


S(nf) 

1 


I fi,i - 1 

K,J + ^i-l,j\ 

Sx 

i,J Ax 

_\ 2 J 

2 / ^1-1,3 J 


when the terms ('^i+i^j + and ^ j + are both positive 

and by 

^(uf) 


3X 


1 

i,J ~ AX 




(30b) 


when the terms 




+ U; 


,,j )/2 “'5 (“ 1,3 “i-1,3)/' 


2 are both negative. 


In Eqs. (30a) and (30b), the function f represents either 9 or co. When 
the average velocities {^>5 

different sign, a mixed expression is required which contains one term 
from each of Eqs. (30a) and (30b), as appropriate. A similar procedure 

is used to approximate ^ and ^^. 7 ^ ^ according to the sign of 

+ v^^jH 2 and + ^i,j-l)/^* product of the average ve- 
locity, for example, + ^i,j)/^ either fj^ j or in Eqs. 

(30a) and (30b) represents the convective transport of f between node 
points (i + 1,0) and (i,j). The selection of f^ ^ or fj.-i according 

to the sign of the mean velocity, is necessary for f to be strictly con- 
served in transport between the nodes. 

Equations (15) and (I6) can be approximated by the following finite- 
difference forms: 


^i,.i ~ ^i,j ^ %^u - 

At Ax Ay 




~ ^^i,j + ^i-l,j 


G 


i»j+l 


^^i,j ^i-l,j 


(Ax)‘ 


(Ay)‘ 


+ 2 


(31) 


and 
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- ‘"i„1 . - “L^L . 

CAt) 


= Pr 


where 


j-^ 


i+l„i ~ •*• ^i-1,3 . ~ ‘^i^j-l 

(A3c)^ (Ay)^ 


- Pr 




i1r - ('^i+i,j ■*■ » 

^L = K,j + > 


% = h,i '"R > 0 

% = ^i+1,3 ^R < ° 

^L = ^i-l»i ^L > 0 

,f“L = ^L < 0 


(32) 

(33) 

(34) 


(35) 



- (n,3+l n,; 

j)/2 

> 

(36) 

Vp = 

(n,j + ^i,j-i)/2 

9 

and 

(37) 

% 

= fi,j ^0^ > 

0 




% 

= fi,J+l 

< 

0 







► . 

(38) 


= f’i,J-l 

> 

0 




= fi,j for Vp < 

0 





It should he noted that the variable f in Eqs. (35) and (38) represents 
either 6 or w. The superscript star in Eqs. (31) and (32) denotes the 
time level n + 1 while all other terms are evaluated at the previous time ■ 
level n. 

By the use of central differences, Eq. (29), to approximate the vor- 
ticity-stream function equation, Eq. (19) > and the velocities in Eq.(l8), 
one obtains 


u 


*’•' ' (Sr)^ ’ '39; 
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9 


and 


(40) 


u 


4 , J - 


_ ^i,3+l ~ 


2W7 


^ - n-l,j 

“ 2l^c) 


(41) 


3 .2 Approximation for the Vorticity at the Boundary 

Equation (32) is used to solve for the vorticity at the interior 
grid points. An explicit relation for the vorticity at rigid boundaries 
is, however, not possible. Therefore, the following procedure was used 
to determine approximate values of vorticity on the boundaries. 

During a time step, the value of toj^j at the boundaries is held con- 
stant. After vorticity values at the interior grid points are calculated 
for the new time step, and the stream function has been computed from it 
at the same time step using Eq. (39) » the values for the vorticity at the 
boundaries are calculated from the stream function. 

The relation for the vorticity at the wall can be derived from a 
Taylor's series expansion of the stream function and the boundary condi- 
tion at the rigid wall. For example, at the lower boundary of the layer, 
\|r(x,o) = 0. It follows then that 


Sx 


y=0 


- ^ 
" Bx^ 


= 0 


y=0 


and from Eq. (l9) j that 


^ y=0 


(42) 


A Taylor's series expansion in y is next used to determine and 

in terms of ti,!’ Noting that (B\lf/By)y_Q, one obtains 

= (43) 




and 




0 ((Ay)^) 

From Eqs. (43) and (44), the value of w at the wall is 


(44) 


( 45 ) 
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Similar relations for the vorticity on the other boundaries can also be 
derived in the same manner. All of the ccmputations done in the present 
study use Eq. (45) and similar equations to evaluate the vorticity at the 
boundaries. An exposition of these equations is presented by Emara [ 4 ]. 

It may be noted that Mayinger et al. , [12], in their study of ther- 
mal convection in an internally heated fluid layer bounded by two con- 
stant temperature boundaries, used an expression of first-order accuracy 
for the vorticity at the walls in the following form: 

^1,1 " (^)g 

They did not use an expression of second-order accixracy, such as Eq. ( 45 ), 
apparently to avoid numerical instability. It is believed that the source 
of the numerical instabilities in the work of Mayinger ^ was the use 
of differencing schemes of different orders of accuracy at the interior 
grid points and at the boundaries. All of the computations done in the 
present study to solve for the vorticity at the wall were obtained by 
using Eq. ( 45 ). 

In order to determine the possible difference between the resvilts 
obtained by using Eqs. ( 45 ) and ( 46 ), a computation was done using Eq. 

( 46 ) . The results of this test computation were in reasonably good agree- 
ment with the second-order method. The difference in the heat transfer 
coefficients at the upper boundary was less than 0.5%. 


3.3 The Computational Method 


New values of the dependent variables over a time step are computed 
in the following manner; 


1. The calculation is begun by setting the initial conditions 

u = v= 0, \|f = u = 0, and 0 = 0 everywhere in the domain of in- 
terest, i.e., 0 g y S 1 and 0 ^ x S X/L. 


2. The time step. At, is determined. This is limited by stability 
considerations which can be briefly explained by the following: 
The differential equations of interest can be written in the 
general form 


3f d^f . df . hf 

at “ a? aF ’ 


('•7) 


where a^ , as, 8 - 3 , and a^ are functions of x, y, and t. An ex- 
plicit finite-difference method to approximate Eq. ( 47 ) reduces 
it to 
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m 


= 


n 


a . -F . 

1,0 


+ a, 




in 


+ a. 




n 




The coefficients aj^j... denote quantities which are constant 
over a time step. The quantities Fj^j... are functions of x, y, 
and t and the superscripts n and n + 1 denote time levels t and 
t + At, respectively. Stability in the sense of Lax and 
Richtmeyer [11] follows if the coefficients ai.j in Eq. (48) are 
positive, i.e., the equations of the "positive" type if 

ai,j ^0 (49) 

for all values of i and j . Another concept defined by Lax and 
Richtmeyer and usually associated with finite-difference equa- 
tions is "consistency." A finite-difference equation [e.g., 

Eq. (48)1, is said to be consistent with the given differential 
equation, Eq. (4?), if. the truncation error involved in replac- 
ing the derivative by finite-differences vanishes as the spatial 
and time increments both approach zero. Lax and Richtmeyer 
showed that if the consistency holds, then stability and conver- 
gence are equivalent, and moreover stability implies convergence 
It is not difficult to show that the present finite- difference 
formulation of the thermal convection problem can be made to 
satisfy the stability requirement, Eq. (49). In Eq. (48), all 
of the coefficients are always positive except aj j, and this 
coefficient can be made positive by restricting t^e size of the 
time step At. The resulting finite-difference equations will of 
the positive type (i.e., stable) provided the following inequal- 
ity holds when Uj^, Uj^, Vy, and Vy are positive: 


At 


|_AX Ay 


2 2 -j-i 

(ax)^ ^ (Ay)^_ 


for the energy equation, Eq. (31), and 


(50) 


for the momentum equation, Eq. (32). Similar restrictions on 
At can be derived when Uj^, Uy, Vy, and Vy are negative [4]. 

It may be noted that if the nonlinear terms in the energy 
equation and the momentum equations are approximated by central 
differences, then the finite-difference equations are ^stable 
when Uj^, Uy, Vy, and Vy are positive, provided that 



2 


+ 


(52) 


and 


At ^ 





for the energy equation, and 


(53) 


At s 



(54) 


and 




< 2 

■" Ax * 




2 Pr 


(55) 


for the momentum equation. Similar restrictions on At can be 
derived when Uj^, ut, Vq and Vp are negative [4]. From the in- 
equalities Eq. (53) and Eq. (55) j one can see the restrictions 
on the grid size which is introduced by using central differences 
to approximate the nonlinear terms in the energy and momentum 
equations. 


3. The temperature distributions at one time step in the future is 
computed using Eq. (31). The results obtained for the tempera- 
ture distribution are then used in Eq. (32) to calculate the 
vorticity at all interior grid points. 


4. Equation (39) is then used to find the stream function at all 

interior grid points. As mentioned earlier, the method of solv- 
ing the vorticity- stream function equation may be different from 
that used in solving the energy and vorticity equations. The 
vorticity- stream function equation is usually solved by an it- 
erative method. The method used in this study is the successive 
over-relaxation method. The iterative relation for this method 
as applied to Eq. (39) is 


4 -i,i = (1 - f2)' 


(ax)^ (Ay)‘ 

n,j+i + ^i,j-l 

(^)= 


>i+i„i +n-i,,i 

_2 L (Ax)^ 


0)1 


, J J > 


(56) 


where Q is the over-rel6ixation factor. The iterations are per- 
formed on Eq. ( 56 ) with 'an algorithm derived from the Gauss- 
Seidel [I 9 ] iteration method. The convergence criterion observed 
in computations is 
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max|\|f?‘*’^ - lif? . 

£jlL £ p 

max|tfl| " ■ 

-*■ 5 «J 

For the present study, = 1.75 and e = 10““^. 

5. The values of the vorticity at the rigid walls are then calcu- 
lated using Eq. (45). 

6. Finally, Eqs. (40) and (4l) are used to compute the velocity 
components u and v at the interior grid points. 

7. Steps 1-7 are repeated until steady state is reached; i.e., the 
difference between the horizontally averaged temperatTire distri- 
butions of two consecutive time steps is of the order of 10“® 
over the entire layer. 


3.4 The Computation of Heat Flux at the Upper Boundary 

The energy flux from the fluid layer to the cooled, isothermal upper 
boundary was obtained from the temperature field. The local Nusselt num- 
ber at the upper boundary was defined as 


Nui = 




The temperature gradient, (57) was obtained by using a 

Taylor's series expansion of the temperature, 0£j, atJ=N, j=N-l, 
and 0 = N - 2. When the terms and are eliminated ftom 

each of these series expansions, the following expression is obtained for 
the temperature gradient at the upper boundary: 


“ 90i,N-l ^®i,N-2 


'i,NW 




From the local heat flux developed in this manner, the average value 
of the Nusselt n\imber was determined by numerical integration using the 
trapezoidal rule. 


L=1 + ^'^l,i=MM f ^ 
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SECTION IV - PRESENTATION AND DISCUSSION OF RESULTS 


4.1 Comparison with Previous Studies 

A comparison with the numerical work of Mayinger et [12] and the 
experimental work of Kulacki and Goldstein [7] was run in order to test 
the precision of the numerical scheme used in the present study. The work 
of Mayinger ^ and Kulacki and Goldstein treated thermal convection 
in an internally heated fl\aid layer confined between two constant-temper- 
ature, rigid plates. In Fig. 2, isotherms and streamlines calculated by 
the method of this study are compared with isotherms and streamlines ob- 
tained by Mayinger ^ al. In Figs- 3-5 » the average temperature distri- 
butions obtained from the numerical scheme of this work are compared to 
the spatially averaged temperature distributions obtained experimentally 
by Kulacki and Goldstein. In both cases, the agreement is good. In Fig. 
6, Nusselt number versus Rayleigh number correlations obtained with the 
n-umerical method of this study and the correlations of Mayinger at al . 
and Kulacki and Goldstein are presented. The agreement between the three 
correlations is seen to be very good, and the numerically derived corre- 
lations of the present study lie, within the scatter of the data in the 
experimental correlation of Kulacki and Goldstein. 


4 .2 Results of the Present Study 

Flow fields, temperature fields, and average heat transfer coeffi- 
cients at the upper surface were obtained for 5 x 10® S Ra g 5 x 10®, 

0.05 g Pr S 20, and 0.125 g L/x g 1 in the present study. The majority 
of the computations were done with layer aspect ratios of 1.0 and 0.5* 

The results presented in this section are for steady convection in ac- 
cordance with the criterion for sl^ady state solutions desc rib e d in i 

: Section 3»3« All of the results were obtained using a grid size 

of Ax = Ay = 0.033j which was chosen based on a consideration of accuracy 

and computation time [4]. 

In Figs. 7-20, typical streamline and isotherm patterns are pre- 
sented for Pr = 6.5 arid several values of Rayleigh number. The results 
for each value of Ra are presented at the same dimensionless time from 
the initiation of the computation to facilitate comparison. It should be 
noted that the results presented in Figs. 7-20 are not truly steady. 

This characteristic vinsteadiness of thermal convection in fluid layers 
with uniform volumetric energy sources has been observed in experiments 
[7,12] and also in the present and previous theoretical studies [12]. 

At a Rayleigh number of 5 x 10'^ (Ra/Ra^ = 36 where Ra = I386 [10]), 
it is seen in Fig. 7 that two counter-rotating swirls (i.e., vortex-like 
flows) occupy the fluid layer. This flow acts to carry warm fluid upward 
over a relatively broad region in the center of the layer (Fig. 8). Cool 
fluid is carried across the upper boundary and down the side walls of the 
layer . It may be noted that the cool down-flows and the warm up-flow 
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T-T,/ (HLV8k) 

Horizontally averaged temperature distribution of Kulacki and 
Goldstein [?] and the present study for a layer with two rigid 
isothermal horizontal boundaries. Ra = 5-13 x lO'^ 





0.2 


T-T,/(HLV8k) 


Fig. 4 - Horizontally averaged temperature distribution of Kulacki and 
Goldstein [?] and the present study for a layer with two rigid 
isothermal horizontal boundaries. Ra = 4.04 x 10® 











Fig. 7 - Streamline pattern.! Ra = 5 x lO'^ and Pr = 6.5 



Fig. 8 - Isotherm pattern.' R^ = 5 x 10“^, Pr = 6.5? and A0 = 0.065 
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Fig. 9 " Streamline pattern. I Ra = 10^ and Pr = 6.5 



Fig. 10 - Isotherm pattern. Ra = 10^, Pr = 6.5» and A0 = 0.06 
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Fig. 11 - Streamline pattern, i Ra = 5 x 10^ and Pr = 6.5 



Fig. 12 - Isotherm pattern* Ra = 5 x 10®, Pr - 6.5» ^^d A0 - 0.04 
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Fig. 13 - Streamline pattern. Ra = 10® and Pr = 6.5 



Fig. l4 - Isotherm pattern » i Ra = 10®, Pr = 6.5> and A0 = 0.035 
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Fig. 15 - Streamline pattern. ' Ra = 5 x 10® and Pr = 6-5 





Fig. 16 - Isotherm pattern*; Ra = 5 x 10®, Pr = 6.5? and A0 = 0.028 
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and Pr 












occupy about equal fractions the horizontal extent of the layer. These 
flow and temperature fields do not persist, however, as the Rayleigh num- 
ber is increased to a value of 10® (Ra/Ea^ =72). At this Rayleigh nvun- 
ber, the flow field is characterized by a reversal in the direction of 
flow at the center of the layer (Fig. 10) and the formation of four 
counter-rotating swirls (Fig. 9)* The regions of down-flow and up-flow 
still occupy about equal fractions of the horizontal extent of the layer, 
and the up-flow is a relatively broad region corresponding to each pair 
of swirls. The region of down-flow, however, becomes thiner at a Rayleigh 
number of 5 x 10® (Ra/Ra^ 360) as shown in Figs. 11 and 12. At 
Ra = 10® (Ra/Rac - 722) , the flow is still driven by four counter-rotating 
swirls (Fig. 13), and the broad region of warm up-flow (Fig. l4) has a 
horizontal extent somewhat larger than Ra = 5 x 10®. In Figs. 15-18, the 
existence of large-scale eddy motion becomes evident. The streamline pat- 
terns for Ra = 5 X 10® and 10^ have lost some of the regularity observed 
in the streamline patterns at lower Rayleigh numbers. Isotherm patterns 
at these same Rayleigh numbers show that thermals released from the ther- 
mal boundary layer at the upper surface have a scale of the order of the 
layer depth. For both of these Rayleigh numbers, down-flows occirr in rel- 
atively narrow regions except when a thermal is released from the upper 
surface (Fig. l8), and up- flows occupy more of the horizontal extent of 
the layer. Thus, at these Rayleigh numbers, down-flows are considerably 
higher in velocity than the up-flow. At a Rayleigh number of 10® (Figs. 

19 and 20), these feattires are even more pronounced. 

In Figs. 21-28, isotherms and streamlines for Pr = 1 at different 
Rayleigh niunbers are presented. The isotherms and streamlines presented 
in Figs. 21-24 at Ra = 10® and 5 x 10® show a form similar to those pre- 
sented in Figs. 9-12 at Pr = 6.5 and the same Rayleigh numbers. At these 
low Rayleigh numbers, the gross features of the flow and temperature 
fields are apparently independent of the Prandtl number. Comparing the 
streamlines and isotherms for a Rayleigh number of 10® and Pr = 1 (Figs. 

24 and 25) with those obtained with Pr = 6.5 and at the same Rayleigh 
number (Figs. 13 and l4) , one can see that although fo\ar swirls exist in 
both cases, a slight difference appears in the feattires of the isotherms. 
In Figs. 26 and 27 are shown the streamlines and isotherms for a Rayleigh 
number of 5 x 10®. These isotherms and the streamlines can be best com- 
pared with those presented in Figs. 15 and I6 for Pr = 6.5* While at a 
Prandtl nimiber of 6.5, three tongues of down-flow exist in the layer; at 
Pr = 1 only two such tongues appear. At the larger Prandtl number, the 
convective motion splits into more swirls than at lower Prandtl numbers. 

To gain further insight on the influence of the Prandtl number. 

Figs. 29 and 30 contain the results for Pr = O.O5 and Ra = 5 x 10®, and 
Figs. 30 and 3I contain the results for the same Rayleigh number but for 
Pr = 20. A Rayleigh ntunber of 5 x 10® has been chosen because at this 
value the difference between the convection form for Pr = 6.5 and Pr = 1 
is significant. At Pr = O.O5, the isotherm pattern (Fig. 30) indicates 
that conduction is a dominant mechanism for energy transport near the up- 
per boundary, more so than at higher Prandtl numbers (compare Fig. 30 
with Fig. 32). At the lower Prandtl numbers, the flow is characterized 
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Fig . 21 - Streamline pattern . j Ra = 10^ and Pr = 1 



Fig. 22 - Isotherm pattern, | Ra = 10^, Pr = 1, and A0 = 0.077 
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Fig. 24 - Isotherm pattern^ Ra = 5 x 10^, Pr = 1, and ^ = 0.04 
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35 





36 


and Pr 









by the absence of narrow regions of down-flow from the upper surface and 
fewer regions (swirls) of circiilating flow (compare Fig. 29 with Fig. 3l) 
A comparison of Fig. 30 with Fig. 32 reveals that the thermal boundary 
layer at the upper surface is much thicker at the lower Prandtl ntunber. 

A comparison of Fig. 29 with Fig. 31 shows, also, that at Pr = 0.05 the 
distance between streamlines in the down-flow zones is about the same in 
the regions of up-flow. This means that in both regions, the velocities 
are about the same. On the other hand, at Pr = 20 it is clear that the 
velocity in the down-flow region is much higher than in the region of 
up-flow. 


The temperature fields for steady convection enable one to calculate 
the average heat fl\uc at the upper boundary and the horizontally averaged 
temperature distribution in the layer. From these quantities, it is pos- 
sible to determine the overall Nusselt number at the upper boundary. It 
was found that the Nusselt n\miber and Rayleigh number could be correlated 


Nui = 0.420 Ra°*^^® 

5 X 10^ S Ra g 5 X 10® 


Pr = 6.5, 0.125 g L/X g ij 


(60) 


and 

Nui = 0.477 
5 X 10® g Ra g 5 X 10® 
0.05 g Pr g 20 


(61) 


0.125 ^ l/x si J 

Agreement between the correlation obtained numerically, Eq. (60), 
and that obtained experimentally by Kulacki and Emara [8] is very good 
(Fig. 33)- The n\mierically derived correlation differs by less than 2^ 
from the experimental results for low Rayleigh numbers (i.e., 

10 g Ra/Rag g 38) and both correlations are nearly identical at higher 
Rayleigh numbers, A comparison between the numerical results for the 
heat transfer coefficient of this study, Eq. (60), and nmerical work of 
Thirlby [16] using the method of artificial compressibility and Roberts' 
analytical results [I5] are presented in Fig. 34 along with the experi- 
mental data of Kulacki and Emara [8]. It can be seen in Fig. 34 that the 
results of the present study give better agreement with the experimental 
data of Kulacki and Emara than either of the : other studies . 


Horizontally averaged temperature distributions across the layer 
obtained numerically in this study are compared with the measurements of 
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Fig. 3^ - Nusselt numbers at the upper boundary from the present study and the studies of Roberts [15 Ij 
Thirlby [16], and Kulacki and Emara [ 9 ] 


Kulacki and Emara [9] in Figs. 35-37- Good agreement is seen between the 
numerical and experimental results. The deviation between the experimen- 
tal and numerical temperature profiles in the core region (Figs. 36 and 
37) is believed to be caused by the distiirbances of the flow introduced 
by a thermocouple probe which was used to obtain time-averaged measure- 
ments of temperature within the layer. Nimierical results for the hori- 
zontally averaged maximum temperature difference across the layer are 
compared in Fig. 38 with the experimental results of Kulacki and Emara. 
The agreement between the measured and computed values is seen to be good 
over the range of Rayleigh numbers considered in the present study. 
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• Experimental 
— Numerical 



Fig. 35 - Horizontally averaged temperatijre distribution determined nu 
merically and the time averaged measTirements of Kulacki and 
Emara [9] for Ra = 5*59 x 10'^ 
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Fig. 37 - Horizontally averaged temperatirre distribution determined riu' 
merically and the time averaged measiorements of Kulacki and 
Emara [ 9 ] for Ra = 8.55 x 10^ 








SECTION V - CONCLUDING REMARKS 


The present study has resulted in a n\jmerical solution of the partial 
differential equations which govern thermal convection with uniform volu- 
metric energy sources in a horizontal fluid layer. The boundary condi- 
tions of primary concern are a rigid, isothermal upper bo-undary and a 
rigid, zero-heat-flvix lower boundary. The side walls are assumed to be 
rigid and perfectly insulating. Subsidiary calctilations have been done 
for a fluid layer with two rigid, isothermal horizontal boxindaries with 
both side walls rigid and perfectly insulating. These calculations serve 
as a check on the numerical formulation of the problem owing to the exist- 
ence of both experimental and ntimerical studies for a layer with these 
boxmdary conditions. 

The numerical results for the temperat\ire distribution within the 
layer permit the computation of the average heat transfer coefficient at 
the upper boundary. A correlation of the Nusselt number in terms of the 
Rayleigh and Prandtl nxmibers has been determined in the following form; 

Nui = 0.477 Ra°-2iOpr°-°^°'^ 

5 X 10^ g Ra £ 5 X 10® 

0.05 ^ Pr S 20 

0.125 g L/X SI. 

When this form of correlation is restricted to the results for Pr = 6 . 5 » 
predicted Nusselt numbers are in excellent agreement with the measure- 
ments of Kulacki and Emara [8]. The resiolts of the present study are thus 
confirmed by the experiments for Ra g 5 x 10® [Eq. (60) and Fig. 33] • 

The results of the present study are formally restricted to two- 
dimensional laminar convection. While the horizontally averaged temper- 
ature profiles and the correlation for the Nusselt number derived from 
them are in agreement with the measurements of Koolacki and ^ara [8,91 > 
it should be noted that the Rayleigh number range considered overlaps what 
is believed to be the turbulent regime of flow. T^e. early experiments of 
Tritton and Zarraga [l8] on the planfo^m of the motion in a layer with 
boundary conditions similar to those of the present study showed that tur- 
bulent motion begins at a Rayleigh nxunber on the order of 80 times the 
critical value for the onset of convection predicted by linear stability 
theory. More recent experiments by Kulacki and Goldstein [7l on a layer 
with two isothermal boundaries indicated that turbulent motion begins at 
Rayleigh numbers about 100 times the critical value. The work of Tritton 
and Zarraga, however, is considered as the only substantiated indication 
of the onset of turbulence in the system of interest here. The critical 
Rayleigh number given by linear stability theory is Ra^ = I 386 [10], and 
this places the results of the present study in the range 3 .61 ^ Ra/Ra,, 

S 3.6 X 10^. Thus the general validity of the isotherm and streamline 
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patterns for high Rayleigh n\xmbers of the present study can be questioned 
despite the agreement between computed and measured Nusselt numbers. It 
is, indeed, possible that overall heat transfer results are relatively 
insensitive to the details of the flow field within the layer as long as 
the scale of the motion (e.g., see Figs. 15-20) is of the order of the 
layer depth. A recent study by Cheung [3]> who developed a semiempirical 
model for turbxilent convection, supports this contention. On the other 
hand, numerically induced viscosity effects may essentially inhibit the 
destabilizing effects of turbulence in the finite-difference computations. 
These and other issues related to the validity of finite-difference cal- 
cxilations of thermal convection over a wide range of Rayleigh numbers can 
be resolved only after further analytical and numerical investigation. 

, In its present form, the computer program developed in this s tudy 
will "enable computations of thermal convection to be carried out 
for a variety of thermal and hydrodynamic boundary conditions of interest 
to both fluid dynamicists and engineers. For example, it is of interest 
to determine the influence of a free upper surface on the heat transfer 
coefficient at the upper boundary and on the overall flow field within 
the layer. Additional computations on the effect of conducting side walls 
or, equivalently, horizontal mean temperature gradients greater than zero 
on the heat transfer rate to the upper boundary would be of interest in 
the fields of geophysics and nuclear safety engineering. These calcula- 
tions are presently being carried out, and the results will be forth- 
coming . 
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